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Abstract — The present work investigates the efficiency of the Multigrid method when applied to solve two- 
dimensional laminar steady free convection flow in a square enclosure partially heated from below. The numerical 
method includes finite volume discretization with upwind scheme on structure orthogonal regular meshes. The 
performance of the correction storage (CS) Multigrid algorithm is compared for different numbers of sweeps in each 
grid level. Up to two grids, for both Multigrid V- and W- cycles, are presented. The results are mainly analyzed in 
terms of the average heat transfer at the walls of the enclosure and Multigrid performance on the rate of 
convergence. It is also shown that convective heat transfer has a characteristic behavior for each boundary 
conditions adopted in given ranges of the governing parameters. 
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Abbreviated title: Numerical Analysis of Free Convection Using the Multigrid Method 

Nomenclature 


C. 

Specific heat at constant pressure 

CPU 

CPU Time (s) 

g 

Gravitational acceleration 

Gr 

Grashof number 

h 

Average convective heat transfer coefficient 

H 

Height of the enclosure 

L 

Domain length / height 

k 

Thermal conductivity of the fluid 

M 

Maximum grid number 

Nu 

Average Nusselt number. 

P 

Thermodynamic pressure 

Pe 

Peclet number 

Pr 

Prandtl number 


Residue 

T 

Temperature 

K 

Temperature of the isothermal vertical wall 
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Source term for^ , (p = U,V, p,t 


u 

Component of velocity along X axis 

V 

Component of velocity along y axis 

w 

Width of the inlet, and the vent 


x,y 

Subscrit 

Cartesian coordinates 


Uj 

Nodal index 


in 

Input values 


k 

Grid level 


nb 

Neighboring 


Greeks Characters 


a 

Thermal diffusivity, pC / k 


P 

Coefficient of thermal expansion. 

.Py 

^dp'' 

S 

Dimensionless length of the heat source. 


V 

Kinetic viscosity of the fluid 



Dynamic viscosity 


p 

Density of fluid 


(p 

General variable 



Diffusion coefficient for^ ,(p = U,V, p,t 


Number of Coarsest-grid iterations 

V Number of pre-smoothing iterations 

V Number of post-smoothing iteration 


I. INTRODUCTION 

This study promotes a discussion surrounding the 
efficiency of the Multigrid method applied in a specific 
configuration. Multigrid methods have been used in many 
different calculations as a result of its facilitated converge 
capability.In single-grids, convergence rates solutions are 
greater in the beginning of calculations, reducing this 
sensibility while the iterative processes goes on. This hard- 
to-converge behavioris due to the iterative methods which 
smooth out only those Fourier error components of 
wavelengths smaller than or equal to the grid size. Naturally, 
this effect becomes more significant as the mesh becomes 
refined. At another level. Multigrid methods cover a broader 
range of wavelengths trough relaxation on more than one 
grid, making it easier to converge. 


A Multigrid cycle is a repetitive procedure practiced at 
each grid level according to the grid hierarchy. The V- and 
W- cycles are types of Multigrid that determines the 
convergence criterion and the number of iterations in each 
step along consecutive grid levels visited by the algorithm. 
Within each cycle, the intermediate solution is relaxed before 
(pre-) and after (post-smoothing) the transportation of values 
to coarser (restriction) or to finer (prolongation) grids [[1]- 

[3]]. 

The Multigrid method can be roughly classified into 
two major categories. The CS formulation, where algebraic 
equations are solver for the corrections of the variables and 
the Full Approximation Storage (FAS) scheme, where 
variables themselves are handled in all grid levels. Since 
much work has been done on both major classifications, 
specific recommendations are admitted. The application of 
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the CS formulation is recommended for the solution of linear 
problems being the FAS formulation more suitable to non¬ 
linear cases [[l]-[3]].Therefore, an exception has been 
disclosed in the work of [[4]], who reported predictions for 
the Navier-Stokes equations successfully applying the 
Multigrid CS formulation. In the literature, however, not too 
many attempts in solving non-linear problems with Multigrid 
linear operators are found. 

Acknowledging the advantages of using multiple grids, 
[[5]] presented numerical computations applying this 
technique to recirculating flows in several geometries of 
engineering interest. There, the correction storage (CS) 
formulation was applied to non-linear problems. Later [[6]], 
analyzed the effect of Peclet number and the use of different 
solution cycles when solving the temperature field within 
flows with a given velocity distribution. In all those cases, 
the advantages in using more than one grid in iterative 
solution was confirmed, furthermore, [[7]], introduced the 
solution of the energy equation in their Multigrid algorithm. 
Temperature distribution was calculated solving the whole 
equation set together with the flow field as well as 
uncoupling the momentum and energy equations. A study on 
optimal relaxation parameters was there reported. More 
recently, [[8]] analyzed the influence of the increase of 
points of the mesh and optimal values of the parameters of 
the Multigrid cycle for different geometries. Also, [[8]-[ll]], 
presented a study on optimal convergence characteristics in 
solution of conductive-convective problems. 

Much work has been done with enclosure geometries, 
studying heat and mass transfer, (simultaneously or not) 
because of its engineering response value. The reason for the 
attention behind the physical nature of buoyance-induced 
flows is well represented in [[12]]. Cooling of electronic 
devices before its excessive heating, recovery of remnant oil 
in petroleum reservoirs, dispersion of atmospheric pollution 
and its implications in adjacent cities,spreading chemical and 
nuclear waste in soil, are just some examples of its 
importance. 

The current work considers that free convection 
conditions can be imposed inside a square cavity and aims to 
study the interactions between buoyancy forces and heat 
elements inside. This application involves the work showed 
in the literature that has been discussed by many authors. The 
interaction between buoyancy forces and the heated elements 
and the numerical analysis of Multigrid solution applied into 
momentum and heat transfer forms the main objective of 
current work. 


II. GOVERNING EQUATIONS 

The following equations emerge from the mathematical 
descriptions of fluid flow and convective heat transfer in the 
enclosure. These governing equations are based on two- 
dimensional, incompressible, laminar flow in Cartesian 
coordinate system. 


dU dV ^ 
dx dy 

dx dy 


1 dp 
p dx 


+ vV"t/ 


( 1 ) 

( 2 ) 


dx dy p dy 


1 dp 


U — + V— = aV^T 
dx dy 


(3) 

(4) 


Where U and V are the velocity components in X and 
y directions respectively, p is the density of the fluid, p 
is the total pressure and V is the kinematic viscosity of the 
fluid. The gravity acceleration is defined by g and /]j. is the 
thermal expansion coefficient. T and are the 

temperature and the reference temperature, respectively, and 
GC is the thermal diffusivity. 


The transport dimensionless parameters, Grashoff ( Gr 
), which provides the relationship between fluid buoyancy 
and viscosity, Prandtl (Pr), that provides the relationship of 
momentum diffusivity and thermal diffusivity and the 
Rayleigh number ( Ra ), which is an associated number of 
buoyancy-driven flow (natural convection) are given by: 


Gr = 




,Pr = 


V 


V 

a 


(5) 


Ra — Gr- Pr 


( 6 ) 


III. NUMERICAL MODEL 

The solution domain consist on a number of rectangular 
control volumes (CV), resulting in a structure orthogonal 
non-uniform mesh. Grid points are located according to a 
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cell-centered scheme and velocities are store in a collocated and internodal distances is sketched in Fig. Error! 
arrangement[[13]]. A typical CV with its main dimensions Reference source not found.. 



Fig. 1- Control Volume for discretization. 


Writing equations (1), (2), (3) and (4) in terms of a 
general variable g) = ^,U,V,T^ with 

one gets, 


A 

dx 


PU(P-T,^ 

V 


d(p'^ 
dx j 



( 

V 


dcp^ 
dy j 



(7) 


After integrating it over the CV of Fig. Error! 
Reference source not found.. 


they are not repeated in this paper. To summarize, convective 
terms are discretized using the upwind differencing scheme 
(UDS) and diffusive fluxes make use of the central 
differencing scheme (CDS). 

The final discretization for grid node P is done by using 
the integrated transport equation (8) with the substitution of 
all approximate expressions for interface values and 
gradients. 

UpCpp = UpCpp -f -f h (9) 

With the east face coefficient, for example, being 
defined as: 

«£ =max[-C^,0]-f (10) 




^{pU(p)F^{pY(p) 


dx 


dy 


dV^j 


SV 


d ( dcp'^ 
’’’ dx 


dx 


V 


d 

+ — 

dy 


In (10), D^ = fx^d^jtxx^ and = (/7(7are 
^theS^hiiive and convective fluxes at the CV east face. 


r — 


dve and 


S 

a^usual, the operator max[a,Z?] 


returns 


( 8 ) 


the greater between a and b . 


A set of algebraic equations are resulted from the 
integration of terms in (8), as one can name it, convection, 
diffusion, and source. Since these procedures are described 
somewhere else (e.g. Error! Reference source not found.) 


IV. MULTIGRID TECHNIQUE 

The wanted effect of the Multigrid approach in this 
work is for the convergence rate to become independent on 
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the grid spacing and the numerical solution faster. The task 
behind the Multigrid method is to involve a hierarchy of 
successively coarsened grids into the iterative solution 
process. When done, an adequate strategy for the movement 
through the different grid levels should follow, along with 
consistently transferring data with the discretization scheme 
between the grids. This process allows an efficient error 
reduction over a wide spectrum of frequencies. 

If an iterative scheme as the one described below is 
applied to the system of equations on a given grid, it turns 
out only those frequencies of the solution error can be 
reduced efficiently, which corresponds to the grid spacing. 
The high frequencies of the error are reduced a few 
iterations, while the low frequencies nearly remain 
unchanged. At another level, the steps usually taken in 
Multigrid algorithm are the reduction of high frequency 
errors(smoothing), computation of residual error(residual 
computation), decimation of the residual error to a coarser 

i L > 


7 < or . Idiabatii 

* > 




Tu 


Adiabatic _ Adiabatic 


, f ► 


grid(restriction), and the interpolation into a finer grid. The 
present development about the Multigrid technique is also 
presented in [[5],[6],[10],[11]] and for this reason the 
development is not repeated here. 


V. RESULTS AND DISCUSSION 

The computer code was run on an IBM PC machine 
with an INTEL CORE 2 DUO 2.0 GHz processor. Grid 
independence studies were conducted such that the solutions 
presented herein are essentially grid independent. Eor both 
cycles, pre- and post-smoothing iterations were 
accomplished via the Gauss-Seidel algorithm while, at the 
coarsest-grid, the TDMA method has been applied [[13]]. 

The Eig. 2 a) represent general geometries that were run 
with the finest grid having 66x66 grid points highlighted in 
Fig. 2 b). 



Fig. 2 - a) Geometries and boundary conditions and b) computational grid. 


The main difference between the present work and the 
work from [[14]] is that here it is used the prescribed values 
(temperature - T^) while in the reference work is used heat 
fluxes. This particularity reinforces a permanent regime with 
constant physical properties of the flux. 


In order to understand better the implications of this 
new configuration and construct a simpler idea of the square 
cavity, for comparison meaning,Eig. 3 shows the streamlines 
and isotherms, respectively, of a clear square cavity heated 
on the left and cooled from the opposing side for Rayleigh 

numbers ranging from 1 x 10^ to 1 x 10® . 
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Steam Functions 


Temperature 


Ra^l* 10^ 


a) 


c) 


e) 


Ra — 
g) 






Fig. 3 - Natural Convection in Square Cavity from left to bottom Ra = 1 * 10^ , Ra = 1 * lO"*", 


1 * 10 ® 

h) 


Ra = 1 * 10^and Ra = 1 * 10®, respectively. 


In Fig. (3a), Ra = 1 x 10^, the streamlines indicates an 
existence ofone centered vortex while corresponding 
isotherms. Fig. (3b), indicates a conductive heat transfer, 
expressed by the almost parallel pattern with the heated wall. 
The vortex is generated by the horizontal temperature 


dT 

gradient across the section. This gradient,-is negative 

dy 

everywhere, inducing a clockwise oriented vorticity. 

When the Rayleigh number is increased to 
Ra =1x10^, Fig. (3c), the vortex in the middle of the 
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cavity starts to shape differently, slightly into a more elliptic configuration is undone, especially in the middle of the 

configuration. The isotherms. Fig (3d), therefore, have a cavity.Temperature gradients are stronger near the vertical 

considerable convection advance, and the parallel walls, but decrease in the center region. 



a) 



b) 




c) d) 


Fig. 4 - Temperatures residues history for different values of the Rayleigh number: a)Ra = 1 * 10^, 
b) Ra — 1 * 10^, c)Ra = 1 * lO^and d)Ra = 1 * 10^ 
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For 7?a=lxl0^, Fig. (3e), the elliptical behavior 
continuous and the centered vortex is stretched so that two 
secondary vortices can be shown inside of it. Heat transfer by 
convection in the viscous boundary layer alters the 
temperature distribution to such an extent that temperature 
gradients in the center of the domain are close to zero. Fig. 
(3e) shows that, with this change in the sign of the source 
term negative, vorticity is induced within the domain. This 
also causes the development of secondary vortices in the 
core. Fig. (3f) continuous its convection advance in the 
square cavity indicating a faster movement of the flux closer 
to the walls. 

Finally, in Fig. (3g), Ra =1x10^, one can see an 
initiative attempt for a three vortices configuration inside the 
main vortex. Corresponding Fig. (3h) shows that heat 
transfer is mostly convective, again due to the faster 
movement on boarders. 


Fig.4 above, shows the residue history for temperature 
with different values of Ra =1x10^ to 5x10^, up to 3 
grids, for the V- and W-cycles. For a three grids, with 
Ra =1x10^ one can notice a slight advantage in using the 

W-cycle, while in = 1 x 10^ that advantage is almost 
unnoticeable. Now, looking at higher Rayleigh numbers, 
Ra =1x10^ to Ra =5x10^, the V-cycle can be better 
stipulated. This change of cycle due to the increase in the 
Rayleigh number happens because of the changing flux 
attempt to turn turbulent. An explanation for this matter are 
on many references of Multigrid approximation on turbulent 
flow, which in most of the cases use the FAS formulation. 

Another configuration, shows the results of 
isotherms and streamlines of a clear square cavity heated on 
the bottom and cooled from the top for a Rayleigh number of 

Ra = 4 X 10^ . The work of [[12]] is used for comparison. 


Reference [[12]] 



>- 



Present Results 



Fig. 5 - Isotherms and Streamlines for a clear square cavity heated from bottom and from the ceiling for 
Ra = 1 * lO"*", comparison between [[12]] and present research. 


In Fig. 5, it is easy to see the similarities in the results of 
bothworks. Studying these results, one can note the plume 


shaped structure starting at the bottom of the cavity moving 
towards the top. Circulatory motion brings the bottom hot 
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temperature stream up to the top wall, substantially 
penetrating into the flow core. 

From these results, is expected a better understanding 
and a certain familiarization with the final configuration. Fig. 
2a, due to the heated bottom in the square cavity. 
Additionally,Fig. 5 representsthe only possible solutions for 
that Ra number, reinforcing the final results and the 
geometry qualitatively. 


Adiabatic 


Tc 


Tc 


Til 


Adiabatic 


Adiabatic 


Final configuration, square cavity partially heated from 
above, is separated into two different forms, denoted as Cj 
and C 2 ■ The Cy structure hasboth lateral walls cooled at 
temperature and top wall adiabatic, while C 2 hasleft 
lateral wall and top wall cooled at temperature and right 
lateral wall adiabatic. Fig. 6 is a representative scheme of 
Cj and C 2 . Notice that for allconfigurations the total 
surfaces of the cooled walls are identical. 


Tc 


Tc 


Adiabatic] 


Th 


Adiabatic 


Adicd^tk^ 


Configuration Cl Configuration C2 

Fig. 6 - Thermal configurations of the cavity. 
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Ra =4x10^^ 



Fig. 7 - Isotherms, vectors and isotherms in the cavity for 6 — 0.03. 


Final results can be seen in Fig. 7, 8 and 9 for Cj 
and Cj conditions, side by side, as Rayleigh number 
increases top to bottom. As S was incremented, the 
streamlines and isotherms related are plotted. 

First observation lies on the boundary thermal 
conditions. When these conditions are symmetrical, about the 
vertical mid plane, i.e., case Cj, the fluid motion is also 
symmetrical and two counter-rotating cells are formed in the 
cavity. The isotherms are also symmetrical about the vertical 
mid plane and it is noticed that the temperature gradient 
becomes steeper and bigger at the hot surface where a 
thermal plume may be located. In case C 2 , fluid flow takes a 
direction towards the adiabatic wall on the side of the cavity. 


making the flow asymmetric and characterized with a single¬ 
cell of anticlockwise circulation. 

In each case the stagnation point is observed at the 
middle of the bottom wall. It is easy to see the similarity 
between case Cj and the bottom-heated arrangement in Fig. 
5, making a plume shaped pattern. Similar behaviors of the 
flow and thermal fields are observed at other Rayleigh 
numbers as an increment for the following Figures. Finally, 
Fig. 10 presents temperatures residues history for different 

values of the Rayleigh number, Ra = 4x10^ and 

Ra = I X 10^ for the case it was presented in Fig. 7. Once 
more again the Multigrid solution has the best results at least 
in the computational effort reduction. 
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Boundary Condition - Cj 


Boundary Condition - C 2 


Ra = lxlO^ 




Ra=lxlO^ 
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Boundary Condition - Cj 


Boundary Condition - C 2 


Ra=lxlO^ 
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a) 

Fig. 10 — Temperatures residues history for different 
Ra = 1 * 10"*^ and b) Ra — 
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